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^ ■ Abstract 

Ch ■ 

We compute the Lyapunov spectrum and the Kolmogorov-Sinai entropy for a 

^ , moving particle placed in a dilute, random array of hard disk or hard sphere 

r~| , scatterers - i.e. the dilute Lorentz gas model. This is carried out in two 

O 

ways: First we use simple kinetic theory arguments to compute the Lyapunov 
spectrum for both two and three dimensional systems. In order to provide 
a method that can easily be generalized to non-uniform systems we then use 
a method based upon extensions of the Lorentz-Boltzmann (LB) equation to 
include variables that characterize the chaotic behavior of the system. The ex- 
tended LB equations depend upon the number of dimensions and on whether 
one is computing positive or negative Lyapunov exponents. In the latter case 
the extended LB equation is closely related to an "anti-Lorentz-Boltzmann 
equation" where the collision operator has the opposite sign from the ordinary 
LB equation. Finally we compare our results with computer simulations of 
Dellago and Posch and find very good agreement. 
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I. INTRODUCTION 



The Lorentz model of a gas of non-interacting particles which collide with fixed scatterers 
has been a basic model for acquiring an understanding of fundamental issues in both the 
kinetic theory of gases and the general theory of non-equilibrium phenomena in fluids and 
solids . In this series of papers we use the Lorentz gas to study important features of the 
chaotic behavior of systems which show typical equilibrium and non-equilibrium behavior 
such as the existence of a spatially homogeneous equilibrium state, and normal diffusion of 
moving particles among the scatterers. Since the particles do not collide with each other 
the Lorentz gas can readily be analyzed in terms of the properties of one moving particle 
in the set of fixed scatterers. Here we consider the scatterers to be placed at random in 
space, subject only to the condition that they are not allowed to overlap with each other. 
The chaotic properties of a particle moving in a periodic array of non-overlapping hard 
disks have been studied extensively, especially for the case where the density of scatterers is 
sufficiently high that the moving particle is unable to travel unimpeded through the lattice 
(the case of finite horizon) . It is known that the random Lorentz gas is a K-system , and 
that the periodic Lorentz gas with finite horizon is a Bernoulli system [|],^. These results 
are sufficient to prove that the gas has a well defined equilibrium state, and that suitably 
defined initial ensemble distributions will approach equilibrium distributions for long enough 
times. However for the random case there are very few analytic results for quantities that 
characterize the chaotic behavior of the moving particle. There is a conjecture by Krylov that 
the positive Lyapunov exponents for the moving particle are proportional to na'^~^v \\i[na'^]~^ 
if na'^ ^ 1 where n = N/V is the number density of scatterers of radius a in volume V, 
and V is the constant speed of the moving particle ||^. This conjecture has been verified by 
Chernov, who argued that for low enough densities the periodic and the random Lorentz 
gas should have the same value of the Kolmogorov-Sinai entropy, and then calculated this 
quantity for a periodic system at low density. Chernov obtained the results ^ 

hxs = A"*" ~ —2nav\n.na^ for d = 2 (1) 
2 



hxs = ~^ '^t — —Svrna^f In na^ for d = 3. (2) 

Here A"*" denotes a positive Lyapunov exponent. Simple considerations of the number of de- 
grees of freedom and the conservation of energy show that for a two dimensional Lorentz gas 
there can be no more than one positive Lyapunov exponent, and for a three dimensional gas 
there can be at most two of them. The quantity hxs is the Kolmogorov-Sinai (KS) entropy, 
which for a closed, isolated ergodic system, such as the one considered by Chernov, is equal 
to the sum of the positive Lyapunov exponents, according to Pesin's theorem Chernov's 
results are only the first terms in the density expansion of hxs and up till the present work 
no further analytic results were known either for the density dependent corrections to these 
results or, for three dimensional systems, the individual Lyapunov exponents. 

We have been able to use familiar methods from the kinetic theory of gases to calculate 
the Lyapunov spectrum and the KS entropy for random Lorentz gases at low densities 
||I0HT3|. We can do this for closed, isolated systems as well as for closed systems in a 
magnetic field, open systems (with THESE TWO CATEGORIES ARE THE SAME AS 
FAR AS I KNOW escape of particles), and for systems where the moving particle is charged 
and subjected to an electric field plus a thermostat which maintains a constant kinetic 
energy even in the presence of the electric field. These two latter cases are of particular 
interest because of their importance for methods that relate dynamical quantities such as 
Lyapunov exponents and KS entropies to transport coefficients, in this case the diffusion 
coefficient for the moving particle [0,0 . This series of papers will be devoted to obtaining 



the low density results for Lyapunov exponents and KS entropies under the various situations 
mentioned above. Extensions of these results to higher densities and to other quantities will 
be presented elsewhere. 

In this paper we shall consider dilute, equilibrium Lorentz gases in two or three dimen- 
sions, consisting, respectively, of randomly placed but non-overlapping, fixed hard disk or 
hard sphere scatterers, and a point particle of mass m and speed v moving among them. 
The collisions with the scatterers are taken to be elastic. In subsequent papers we will gen- 
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eralize this model to non-equilibrium situations with thermostatted electric fields, and/or 
open systems with absorbing boundaries. The plan of this paper is as follows: In Section II 
we present an elementary kinetic theory argument which correctly provides the low density 
values of the Lyapunov exponents and KS entropy of the Lorentz gas in two and three dimen- 
sions for closed systems and without external fields. In Section III we consider a more formal 
approach to these quantities due to Sinai P], which will form the basis of the extension of 
the kinetic theory approach to non-uniform systems. There we provide the fundamental 
geometric formulas of Sinai which relate the Lyapunov exponents and Kolmogorov-Sinai 
(KS) entropy of a Lorentz gas to the properties of a radius of curvature matrix. Then, we 
use the ergodic properties of the moving particle to express the Lyapunov exponents and 
KS entropy in terms of averages over an equilibrium ensemble. In Section IV we show that 
the pertinent distribution functions can be obtained from the solution of an appropriate 
extended Lorentz Boltzmann (LB) equation, and in Section V we calculate the KS entropy 
of the two and three dimensional Lorentz gas at low densities. In Section VI we consider 
the negative Lyapunov exponents and show how they can be obtained from a solution of 
an "anti-Lorentz-Boltzmann" equation. This will be important for the extension of the 
theory to treat thermostatted systems. We conclude in Section VII with a comparison of 



these results with the results of computer simulations by Posch and Dellago |I6|, and with 
a discussion of the applications of the methods developed here to more general systems. 



II. VELOCITY DEVIATION METHOD FOR LYAPUNOV EXPONENTS AND KS 

ENTROPIES 

We consider a system of N d-dimensional hard sphere scatterers placed randomly in space 
in a dimensional volume V at low density. Here d = 2,3, the spheres have a radius a, and 
the number density of the spheres, n = N/V satisfies na'^ ^ 1. The moving particle travels 
freely between elastic collisions with the scatterers. The phase point x of the particle, i.e. 
its position and velocity, x = {r,v) satisfies the equations of motion f = v; v = between 
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collisions. At a collision of the moving particle with a scatterer, the velocity of the moving 
particle changes according to 

= V — 2{v ■ n)n, (3) 
where v'^ is the velocity after collision, and n is the unit vector in the direction from the 



center of the scatterer to the point of impact of the moving particle at collision [|l7l- This 
completely specifies the dynamics of the particle given its initial phase and the locations of 
all of the scatterers. 

The Lyapunov exponents characterize the rate of exponential separation or of exponential 
convergence of infinitesimally nearby trajectories on the 2d—l dimensional constant energy 
surface S 0. Since the Lorentz gas is a symplectic, Hamiltonian system, if there are non-zero 



Lyapunov exponents they come in pairs of positive and negative values, ±Aj |T^. However, 
the Lyapunov exponent for displacements in the direction of the trajectory is zero since two 
phase points on the same physical trajectory will follow one another without exponential 
separation or contraction. Therefore there can be at most d — 1 positive and d — 1 negative 
Lyapunov exponents for our system. 

To treat the Lyapunov exponents we consider a bundle of infinitesimally nearby trajec- 
tories on S and follow the motion of this bundle in time. If the phase point on one reference 
trajectory in this bundle is given by x{t) we denote the deviation of another trajectory in the 
bundle from x{t) by 6x{t) = {Sf{t), 6v{t)). Equations of motion for 6x{t) follow immediately 
from the equations of motion for x{t). Since the moving particle has only kinetic energy the 
requirement that both trajectories lie on the same energy surface immediately leads to the 
property v(t) ■ 6v(t) = for all t. Without loss of generality, we may replace 6r(t) by the 
vector of closest approach between the two trajectories, i.e. we set 6f{t) ■ v{t) = (from 
here on we will use the notation 6f{t) for this vector of closest approach of the perturbed 
trajectory to r{t). Notice that if 6f{t) ■ v{t) = at t = it remains so at all later times, by 
virtue of Eqs. (3-6).). In between collisions the spatial and velocity deviations change with 
time according to 



5r = 5v (4) 
6ij = 0. (5) 

The change of 6x{t) at colhsions requires some analysis which has been provided by 
Gaspard and Dorfman and by Dellago, Posch and Hoover ||2^. These authors have 
shown that the change in Sx at the instant of a coUision is given by 

Sf^ = [1- 2hn] ■ 6f (6) 

^ 2 ^„ _ ^2 

SiP' = [1 — 2hn] ■ 6v H — \vn — nv + — -nn — (v ■ n)l] ■ Sr. (7) 

a [V ■ n) 

where v, Sf, 6v are the velocity of the moving particle, the spatial deviation and the velocity 
deviation of the nearby trajectory, immediately before the collision with the scatterer, while 
the "+" variables denote the values immediately after collision. It is important to note 
here that the velocity deviation, 6v, does not change between colhsions but does undergo 
an instantaneous change at each collision with a scatterer. 

Suppose that we prepare a trajectory bundle with initial values of 6x{0) and we follow 
the motion 6x{t) in time. We can relate the largest positive Lyapunov exponent to the 
asymptotic growth of the ratio |(5'i7(t)|/|(5'i7(0)|, by 



Amax = lim 7 In 



\6v{t)\ 



(8) 



\6vm. 

This follows from the observation that if there are positive Lyapunov exponents, then two 
infinitesimally close trajectories will eventually separate in time unless they are so precisely 
arranged that they approach each other exponentially in time. However the latter situation 
occurs only on sets of zero measure (the stable manifold) in tangent space. Furthermore, 
this exponential separation occurs both in configuration space and in velocity space with 
the same exponential factor, since velocities and positions are related by a simple time 
derivative which does not affect the exponential separation rate, and we may consider the 
separation in velocity space alone. In the next section we consider another calculation of 
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the Lyapunov exponents, which treats the separation in configuration space, and we obtain 
the same results. 

Similarly, the sum of all of the positive Lyapunov exponents can be obtained by following 
the growth of a volume element in velocity space as 

^ Ai = lim -ln|detA(t)|, (9) 

with A{t) describing the linear relationship between 5v{t) and ^■^(O), i.e. 

Sv{t)^A{t)-Sv{0), (10) 

This result follows because the time evolution of the vector 5v, which has d — 1 independent 
components, is dominated by the {d— 1) largest eigenvalues and corresponding eigenvectors 
of the matrix A{t), which are precisely the positive eigenvalues. Suppose now that the 
moving particle undergoes a series of Af collisions in the time interval (0, t) with scatterers 
which we label Si, S2, sj^. Since the width of the trajectory bundle is infinitesimal, each 
trajectory within it has the same number of collisions with each scatterer in the same time. 
Since Sv, as noted before, changes only at coUisions, one has 

\6m _ \dv^\ \6v^_,\ \Svt\ 

where Sv^' is the velocity deviation immediately after the collision with scatterer Sj. For the 
same reason A{t) can change with time only at the instants of the collisions of the moving 
particle with the scatterer, so that 

5v{t) = 5v)^ = Bat • Sv^-i = SL^f ■ bat-i • • • ai • Sv{Q). (12) 

Here is a matrix, to be defined below for the case where the density of scatterers is low, 
that expresses the change in the velocity deviation when the moving particle collides with 
scatterer Sj. Consequently, 

Af 

det A(i) = Jldeta^. (13) 



Of course the number of collisions, TV in time t will depend upon t and the initial value of 
the phase point x(0). Expressions for the largest Lyapunov exponent and for the sum of 
positive Lyapunov exponents can be obtained from Eqs. M - O) as 



M I \5vt\ ^ 

Xmax = hm —77 2^ In I and (14) 

t^oo t Al ^ l^^ll 

Af 1 

^ Aj = lim ——:^ln| detail (15) 

To proceed further we need to use the fact that the density of the scatterers is small, 
that is, na'^ <^ 1. This will allow us to determine the mean of the quantities appearing in 
the sums in Eqs. |15|) . Referring to Eq. (|^ we note that 5r appearing on the right hand 



side is the spatial deviation of the moving particle just before a collision with a scatterer. 
Let us suppose that we consider the collision with scatterer Sj. Then immediately before this 
collision Sr^ = 5r^i + Ti6v^_i, where 6xf_i denotes the spatial and velocity deviations just 
after the collision with the previous scatterer, and is the time between the collision 
with scatterer Sj_i and scatterer Sj. At low scatterer densities the time between collisions 
will typically be inversely proportional to the density of scatterers, so that the ratio of the 
order of magnitudes of the first to the second term in the above expression for Sf^ will 
approach zero as the scatterer density decreases. Therefore to leading order in the density, 
6f~ = Ti6v^_i. We then obtain a low density value for given by 

Svt = (1 - 2nn) + ^ (^vt-i^ - nv^i - (^1 • n)l + (^Y'.^) ^^) ■ Sv^i ^^^^ 

Now we have expressions for the change in the velocity deviation at collision and for the 
matrix a, both of which are needed for the calculations outlined above. To evaluate the 



sums appearing in Eqs. (0, [T^), we note that at low densities none of the collisions are 
correlated with any previous collision, that is, the leading contribution to the Lyapunov 
exponents come from collision sequences where the moving particle does not encounter the 
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same scatterer more than once in the sequence Q. Therefore we can treat each term in 
the sums in Eqs. (0, |1^) as being independent of the other terms in the sum. We have 
expressed Xmax and the sum of the positive Lyapunov exponents as arithmetic averages, but 
for long times and with independently distributed terms in the average, we can replace the 
arithmetic averages by ensemble averages over a suitable equilibrium ensemble. That is 



A. 

and 





"|5t7+|" 




_\6v-\_ 



(17) 



^ Ai = i/(ln|deta|) (18) 

Ai>0 

where z/ is the (low density) value of the collision frequency, J\f/t as t becomes large, and the 
angular brackets denote an equilibrium average. 

We now consider a typical collision of the moving particle with one of the scatterers. The 
free time between one collision and the next is sampled from the normalized equilibrium 



distribution of free times ||T^, P{t) given at low densities by 

P(r) = ue^'^ (19) 

The construction of the matrix a requires some geometry and depends on the number of 
dimensions of the system. In any case we take the velocity vector before collision, v to be 
directed along the z-axis, and take n ■ v = —vcos(f), where — 7r/2 < < 7r/2. The velocity 
deviation before collision 6v^ is perpendicular to the 2;-axis. Then it is a simple matter 
to compute and | deta|. For two-dimensional systems 6v and the matrix a are 

given in this representation by 



^In two dimensions the particle will hit the same scatterer an infinite number of times. However 
the effects of such processes are of higher density, and can be neglected here since the times between 
successive collisions with the same scatterer become typically very large as the density of scatterers 
approaches zero. 

9 



5v 



\5v- 



a 



^^(1 + A) cos 20 sin 20 ^ 



(20) 



(1 + A) sin 20 -cos 20^ 
where we introduced A = (2fr)/ (acos0). To leading order in vt /a we find that 

\5v+ 



\5v- 



A; 



det a| = A. 



(21) 



For three dimensional systems the unit vector n can be represented as n = — cos z + 
sin cos a x + sin sin a y. Now the ranges of the angles and a are < < 7r/2 and 
< a < 27r. There is an additional angle ■0 in the x, y plane such that the velocity deviation 
before collision 5v~ = \ 6v~ \ [x cos ip + y sin -0] . It is somewhat more convenient to use a 
symmetric matrix, a = (1 — 2nn) ■ a, given by 



1 + A(cos^ + sin^ cos^ a) A sin^ cos a sin a 
A sin^ cos a sin a 1 + A(cos^ + sin^ 0sin^ a) 
1 



(22) 



One easily finds 



2tv 



\5v' 



— h sm [a — ip) cos 



cos' 



(23) 



and 



1 ~ 1 f'2vT\^ 
det a = det a = ( I 



(24) 



to leading order in vr/a. 

To complete the calculation we must evaluate the averages appearing in Eqs.(|l3, |18|). 
That is we average over the distribution of free times and over the rate at which scattering 
events are taking place with the various scattering angles. Additionally in 3 dimensions an 
average over a stationary distribution of angles ip has to be performed in general. Due to 
the isotropy of the scattering geometry can here be absorbed in a redefinition a' = a — ijj 
of the azimuthal angle a. This will not be true any more if the isotropy of velocity space is 
broken (e.g. by an external field). The appropriate average of a quantity F takes the simple 
form 
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{F) = ^ dr J dn cos (j)P{T)F, (25) 

where -P(t) is the free time distribution given by Eq. (^9|) and J is a normahzation factor 
obtained by setting F = 1 in the numerator. The integration over the unit vector n, i.e., 
over the appropriate sohd angle, ranges over —njl < (p < 7r/2 in two dimensions and 
over < < 7r/2 and < a < 27r in three dimensions. After carrying out the required 
integrations we find that 

A+ = Xmax = 2nav[- ln(2na^) + !-€] + ■■■ (26) 

for two dimensions. Here C is Euler's constant, and the terms not given exphcitly in Eq. 
(p6|) are higher order in the density. Similarly, for the three dimensional Lorentz gas we 
obtain 

A+„, = na'vrri- ln(n/2) + In 2 - ^ - C] + ■ ■ ■ , (27) 

A+ax + A+,„ = 2na'v7r[- ln(n/2) - C] + ■ ■ ■ (28) 
from which it follows that 

A+,„ = na\7r[- ln(n/2) - In 2 + ^ - C] + ■ ■ • , (29) 
where fi = na^n. We have therefore determined the Lyapunov spectrum for the equilibrium 



Lorentz gas at low densities in both two and three dimensions [|I0|,|T2[. We note that the 
two positive Lyapunov exponents for three dimensions differ slightly, and that we were able 
to get individual values because we could calculate the largest exponent and the sum of 
the two exponents. We could not determine all of the Lyapunov exponents for a (i > 3 
dimensional Lorentz gas this way. Moreover, for a spatially inhomogeneous system, such 
as those considered in the application of escape-rate methods, the simple kinetic arguments 
used here are not sufficient and Boltzmann-type methods are essential for the determination 
of the Lyapunov exponents and KS entropies. 
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We will comment further on these results in Section VII after we have obtained them 
again by a more formal method based upon the radius of curvature matrix method of Sinai 

i- 



III. THE RADIUS OF CURVATURE MATRIX 

Our analysis in this Section is based upon the geometric arguments given by Sinai 
for the relationships between the Lyapunov exponents, the KS entropy, and the radius of 
curvature matrix, which describes the time evolution of the separation of nearby phase space 
trajectories of the moving particle. Unlike the method presented in the previous section this 
method treats the time evolution of the spatial separation of a bundle of trajectories rather 
than the evolution of the velocity separation of the bundle. Here we summarize these 
considerations, referring the reader to the literature for further details P,p!9|. 

As before, the trajectory of the moving particle is specified by the phase x{t) = 
(r(t), v{t)), and we take a. 2d — 2 dimensional plane S, (5rj_(t), Sv±{t)), through x{t), where 
both 6f±{t) and Sv±{t) are perpendicular to the velocity v{t). The nearby trajectories will 
intersect S, and we measure the separation of the trajectories by vectors of dimension 2{d—l) 
in E, Sx±(t), given as 



5x±{t) 



(30) 



The time development of 5x±_ is given in terms of a monodromy matrix M(t,to) satisfying 

fex(t) = M(t, to)- 5x^(^0 ), (31) 

The matrix M(t,to) follows the motion of the particle. It changes continuously with time t 
in the intervals between collisions, and undergoes a discontinuous change at the instants of 
collisions of the particle with the scatterers. Between collisions the monodromy matrix has 
the form 
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M{t,t 



J free flight 



( 1 {t-to)l\ 



yO 1 J 



(32) 



At the instant of a collision there is a discontinuous rotation of the velocity of the moving 
particle from its value before collision, v, to its value after collision, = v — 2n{v ■ n), as 
before. Since the velocity of the particle changes discontinuously at collision, the plane S 
also rotates, and the components of the displacement vector 5x^ change instantaneously. 
The changes in the components of 5f±_ and 5v± at the instant of collision are given in Eqs. 
(^, 0) above. 

In order to determine the Lyapunov exponents for this system we need to examine the 
rate of separation or of approach of infinitesimally close trajectories. This can now be done 
with the aid of the radius of curvature operators p„ and p^, acting on a c? — 1 dimensional 
space of velocity deviation vectors orthogonal to v. Here the subscripts u and s denote 
operators describing unstable, or expanding, and stable, or contracting, trajectory bundles 
respectively. The operator p„ is defined by the relation that 

5r^{t) = -pSt)-5v^{t), (33) 

V 

together with the conditions ■ v = v ■ = 0. This definition is motivated by the 
observation that if the velocity deviation 6v describes separating trajectories, then we can 
apply ray optics to describe the separation of the trajectories 0. In the transverse plane 
this separation will be given by an arc length equal to a radius of curvature multiplied by 
an infinitesimal initial angular separation. We have chosen the units in Eq. (^) so that 
the radius of curvature operator has the dimension of a length. The radius of curvature 
operator can be represented as a (rf — 1) x (rf — 1) matrix since both 6'r± and 6v± are defined 
in a plane perpendicular to v. Now suppose that we consider some initial velocity deviation 
6v±{0) corresponding to a diverging pencil of rays, and we want to obtain an equation of 
motion for the radius of curvature matrix p^. We use the fact that the motion of the 
particle consists of periods of free flight punctuated by instantaneous collisions with the 
fixed scatterers. We first consider the free flight motion. From the fact that in free flight 
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Sr±{t) = (5rj_(0) + tSv j_{0), and Sv±(t) = Sv±{0), we infer that during free flight over a time 
interval t from some initial time t = 0, with initial value p(0), the radius of curvature matrix 
changes with time according to 

Puit) = pM + vtl^. (34) 

where 1^ = 1—vv with v a unit vector in the direction of v. Next we use Eqs.(P,|^) to obtain 
the relation between the radius of curvature operator immediately before a collision, p^"-*, 
and its value immediately after a collision, p[^^. This calculation, while straight forward, 
requires a careful analysis in order to obtain a correct expression for p as a (rf — 1) x [d — 1) 
matrix. This analysis is presented in Appendix A. There we find that at a collision p changes 
according to 



[p.^^^^] = u|[p-(-)i-- 



a 



vh + hv — -nn — (v ■ n)l 

[V ■ n) 



U, (35) 



where U is the reflection operator 1 — 2nn. The inverse radius of curvature tensors [p^ ^'•^-'] 
and [pu ^*-^''] are defined on the subspaces orthogonal to v respectively v' and are extended to 
tensors on full space by requiring that their left and right inner products with v respectively 
v' are zero again. 

To obtain the connection between the radius of curvature operator and the Lyapunov 
exponents for the Lorentz gas we proceed as follows. Suppose we wish to know the values 
of the spatial separation of a bundle of trajectories after a sequence of n collisions labeled 
1,2, ...,n that take place at times ti,t2, ...,tn. We use the fact that just before the collision 
at time tn, the spatial deviation vector is given by 

(36) 

= 1± + t^r„,„_ip-i(+)(t„_i)] ■ 6ff\tn-i), 

where Tj_j_i = ti — tj-i. We also use Eq.(|^), namely that 5rj^^ = U ■ 6r]_ \ Then by iterating 
this equation we obtain 
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1± + fT-„,„_lP„^(+)(t„-l) 
1_L + VTn-l,n-2PZ^''^Ktn- 



■U(n- 1)- 
■ U(n - 2) • • • 



(37) 



where the initial time has been set to t = 0, with initial values indicated for the spatial 
deviation vector, Sfj_{0) and for the radius of curvature operator p„(0). Also U(j) = 1 — 
'2'hj'hj is the reflection operator at the j-th collision. It is important to note that det U(j) = 
-1. 

If there is an exponential separation of trajectories, then we would expect that | |5rj_(t) 1 1 ^ 
(exp At) ||5rj_(0)|| for very large t. A volume element constructed at time t = with d ~ 1 
linearly independent vectors 6fj_{0y is also expanding exponentially with time. The ex- 
ponential growth factor is the sum of the positive Lyapunov exponents and satisfies as 
consequence of Eq. 



lim(__,oo T In det 



hm^^oo I Eo"' In det 1^ + t;ri+i,,p-i(+) {t 



(38) 



lim- 



where in the last expression the logarithm of a determinant is expressed as an integral of a 
trace of an inverse matrix. With the aid of Eq. (^) we can express the sum of the positive 
Lyapunov exponents as 



E A.= hm^ r drTTpZ'ir). 
K>o t Jo 



(39) 



Eq. (^) is Sinai's formula for the KS entropy for a moving particle in a system of fixed 
hard-sphere scatterers By combining Eqs. (Q) and (|35|) one may obtain a continued 

fraction representation for [p^^(t)] which, for a fixed final phase point x{t) and initial 
p„(0), converges rapidly with increasing t. So far we have not used any properties of the 
arrangement of the scatterers, so this formula is still quite general. In the case that the 
system is ergodic the time average can be replaced by an ensemble average, taken with 
an appropriate ensemble distribution function, so we can express the sum of the positive 
Lyapunov exponents as 
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J2X, = v{Ttp~'}, (40) 

Ai>0 

where the angular brackets denote an average over an appropriate stationary ensemble 
reached in the course of time from smooth initial distributions. In the case of interest 
here, this distribution will be an equilibrium distribution but in subsequent papers we will 
need to consider more general steady state distribution functions. 

Before completing this section we wish to give a simple derivation of Eq. ([39| ) which 
applies to a Lorentz gas with any reasonable interaction between the moving particle and 
the scatterers, and which will be used often in the subsequent papers. We use the fact that 
having defined the radius of curvature matrix p^, we may write 

= 5Mt) 

= v[pZ\t)]-Sr\it), (41) 

with solution 

Sf^{t) = Texpv [' drlp^ir)]-^ ■ 6f±{0). (42) 
Jo 

Here T denotes a time-ordering operator. Using the method of differential forms or equiva- 
lent methods we see that the growth of a volume element, 6Vr{t), in configuration space 
is given by 

f|^exp(„{..T.|p„Ml-). (43) 

This result leads immediately to Eq.(^) for the sum of Lyapunov exponents. Vattay has 
shown how to construct the inverse of the radius of curvature matrix for a general poten- 



tial p2[. It is straightforward generalizing the results obtained here to other short range 
interaction potentials and it may well be possible to treat some cases with long range inter- 
actions between the scatterers and the light particle. For the case of hard disks or spheres 
considered here, we can use the fact that the radius of curvature matrices at different times 
commute with each other if the times involved are all within the same time interval between 
one collision of the moving particle and a scatterer and the next collision, to write 
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(44) 



j=0 



where 



R(tj+i, tj) = exp V 



■tj+i 



1 



(45) 



and the prime on the product denotes that the times are to be ordered so that the times 
decrease from left to right in the product. By using Eq. ( ^4|) and carrying out the required 



We will use these expressions in Section VI and in the appendices. 

In the next section we discuss the distribution function appearing in the above ensemble 
average for the case that the scatterers are distributed at random with very low density, 
i.e., the mean free path of the moving particle is very large compared to the radius of 
the scatterers. In order to obtain the individual Lyapunov exponents we have to find the 
eigenvalues of the operator that appears on the right hand side of Eq. (p7|). This operator 
can be expressed as a product of {d — 1) x {d — 1) matrices, which describe the collisions of 
the moving particle with the scatterers, and the free motion in between collisions. Again, if 
the system is at low density the product of matrices can be considered to be a product of 
randomly distributed matrices, since the time between collisions and the collision parameters 
will be sampled from a random distribution, corresponding to the random placement of 
scatterers. In Appendix B we show how methods from the theory for eigenvalues of products 
of random matrices ||2^ can be used to obtain the largest Lyapunov exponent. At higher 
densities of scatterers correlations between collisions will arise, and one will have to take 
into account these correlations when computing the eigenvalues of the product of matrices. 



integrals, one can easily see that this expression is equivalent to Eq. (|37|) . Moreover one can 



express the sum of the positive Lyapunov exponents as 




(46) 
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IV. THE EXTENDED LORENTZ-BOLTZMANN EQUATION 



In order to evaluate the ensemble average appearing on the right hand side of Eq. (^OD 
we need to construct an equilibrium, or more generally, a steady state distribution function 
for the radius of curvature matrix. The physics of the problem suggests that a method 
based upon the Lorentz-Boltzmann equation is appropriate here. That is, we have a particle 
moving in a random array of scatterers making only binary collisions with the scatterers. The 
particle moves freely between collisions and at a collision both the velocity and the radius 
of curvature matrix change instantaneously. Methods are now well known P|,|^ JT7|P5|J2B|] 
for obtaining a generalized Lorentz-Boltzmann equation for the time dependent space and 
velocity distribution function, /(r, -u, t), for the moving particle as a function of the density 
of the scatterers, at least in the case that the scatterers are non-overlapping 0. In this case 
it is possible to obtain an equation for the moving particle which includes the effects of 
uncorrelated binary collisions of the particle with the scatterers, excluded volume effects 
and the effects of correlated collision sequences of the moving particle with the scatterers. 
To lowest order in the density of the scatterers, the distribution function f{f,v,t), satisfies 
the Lorentz-Boltzmann (LB) equation |ll|fl JT7|j26[ 



dt df dv 



na J dh\v ■ n\[Q{v ■ h)f{f,v — 2{v ■ n)n,t) — Q{—v ■ h)f{f,v,t)]. (47) 

Here n is a unit vector in the direction from the center of a scatterer to the point of impact at 
a collision, and G(x) denotes the unit step function. The right hand side of the LB equation 
describes the change in / due to collisions as the difference between the gain and the loss 
of particles with velocity v from "restituting" and "direct" collisions, respectively. Higher 



^The case of overlapping scatterers is complicated by the fact that regions may exist where the 
particle would be trapped for all time. In a transport problem these regions need to be treated 
carefully, since particles trapped in them will not diffuse beyond the borders of the trap. 
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order density corrections to the right hand side of Eq. (^T]), comprising the generahzed LB 
equation, can be obtained using the appropriate set of BBGKY hierarchy equations and 



cluster expansion methods p5| , p6| , |24| . These density corrections have been studied in some 
detail, and it is well known that the so-called ring collision sequences are responsible for 
both logarithmic terms in the density expansion of the diffusion coefficient of the moving 
particle and the long time tails in the velocity auto-correlation function of the moving 
particle |27 . 



Our purpose here is to extend the LB equation by including the radius of curvature 
matrix among the variables described by the distribution function for the moving particle. 
We will consider only the low density version of the kinetic theory for this distribution 
function and leave the discussion of higher order density corrections, including the effects 
of correlated collision sequences on the Lyapunov exponents, to later publications. Thus we 
wish to determine an equation, valid at low densities, for an extended distribution function 
F(r, V, p, t), where we dropped the subscript u on p, and we relate the distribution functions 
F and / by 

f{f,v,t) = J dpF{r,v,p,t). (48) 

Given the stationary solution of the extended LB equation for F, we can determine the sum 
of the Lyapunov exponents as 

\ = v f dfdvdp [Trp~i] F(f, v, p), (49) 

Ai>0 

assuming that F is properly normalized. 

An extended LB equation for F that reduces to the usual LB equation for / upon 
integration over the radius of curvature matrix elements, can be obtained by following the 
heuristic derivation of the LB equation and simply modifying it to include the additional 
variables. That is, we consider a large collection of moving particles in the random array of 
scatterers and ask for an equation for the probability that a moving particle has its values 
for f,v,p in the range df,dv,dp about f,v,p, all at time t, i.e., F{r,v^ p^t)dfdvdp. This 
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probability changes in time due to free motion of the particles and due to collisions. The 
change in F due to the free motion of the particles in time dt is 

[F(r + vdt, V, p + vdtl±, t + dt) — F{r, v, p, t)]dfdv dp = 

(50) 

§ + v-§ + vEtl§^] dfdvdpdt. 
We used Eq. (^) to treat the change in the radius of curvature matrix during free parti- 
cle motion, and we have assumed that there are no external forces acting on the system. 
Otherwise we would need to include terms accounting for the changes in velocity and in the 
radius of curvature matrix over a time interval dt due to the external force. If there were 
no collisions taking place in the system, then the right hand side of Eq. (^) would be zero. 
However the collisions account for the fact that the number of particles at f+vdt, v, p+vdtl± 
at time t + dt is not equal to the number of particles at r, v, p at time t. To account for 
the change in F due to collisions we consider the restituting and direct collisions separately. 
The direct collisions result in a loss of the particles with r, v, p over the interval dt due to 
collisions with scatterers. Elementary kinetic theory considerations ]r7|p5| show that this 
loss is 



J dh\v ■ h\Q{—v ■ n)F{f,v, p,t)dfdv dpdt. (51) 



The restituting, or gain term is found by considering those collisions taking place in the 
time interval t to t + dt that produce particles with r, v, p. Again, elementary kinetic theory 
considerations show that this gain is given by 

dp' 



na'^ ^ I dn\v ■ h\ Q{v ■ fi) 
where w ■ n > 0, and 



dp 



F{v',f,p',t)dfdvdpdt, (52) 



if = V — 2{v ■ h)h, (53) 

and the radius of curvature matrix p' is a restituting value such that the radius of curvature 
matrix becomes p after collision. From Eq. ( P^D and the identity U ■ -0 = -0' one obtains the 
relationship 
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—vh — hv + 



[V ■ n) 



-hn + {v ■ h)l 



U 



(54) 



with v-n > 0. For two dimensional systems the radius of curvature matrix can be represented 
by a scalar (namely its non- vanishing eigenvalue), and the restituting value, p', is given by 

(55) 



p p a cos (p 

where is the angle between n and v with — 7r/2 < < 7r/2. For three dimensional systems 
the radius of curvature matrix can be represented by a 2 x 2 matrix by choosing the principal 
axes of the coordinate frame orthogonal to v. Defining angles through cos (p = n ■ v with 
< < 7r/2 and a as the angle between the second coordinate axis and the plane through 
n and v, and multiplying Eq. (Q) from the left and the right by U one can rewrite this 
equation as 



Up'^U = p 



-1 



2 cos( 



^ ^ 








(56) 



with 



^ 1 + tan^ 6 cos^ a tan^ 6 sin a cos a ^ 



tan^ f/)sina cosa l + tan^0sin^a 



(57) 



Here cos0 = n ■ v with < < 7r/2 and a is an azimuthal angle for h in the plane 
perpendicular to v with < a < 2tt. We can now simplify the restituting term in the 
extended LB equation by noting that we can combine the Jacobian with the distribution 
function F to write 



dp' 



dp 



F{p')= I dpi5{p-p{pi)) F{pf), 



(58) 



where the S function in the integrand selects the right restituting value of the radius of 
curvature matrix in accordance with Eq. (^^. Putting everything together, we can obtain 
the extended LB equation as 
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(59) 



dt df ^ dpii 

na'^~^ j dh\v ■ n\ Q{v ■ h) J d p'5 {p - p{p')) F{f,v' , p' ,t) - Q{-v ■ h)F{f,v, p,t) 

In the next section we will use Eq. ( ^9|) to compute the KS entropy for two and three 
dimensional Lorentz gases in equilibrium. Before turning to this calculation, we make an 
observation about the restituting radius of curvature matrix. We note that the diagonal 
elements of the radius of curvature matrix will grow with time between collisions. Thus the 
average value of the diagonal elements of the radius of curvature matrix before collision will 
be on the order of the mean free path between collisions, /. For low density of scatterers, the 
mean free path, /, will be much larger than the radius of the individual scatterers, a, such 
that a// ~ na'^ <^ 1. This observation will allow us to greatly simplify the delta function in 
the restituting collision term in Eq. (|59D, and thereby simplify the calculations to follow. 



V. EQUILIBRIUM SOLUTIONS OF THE EXTENDED LB EQUATION 

In this section we construct the equilibrium solutions of the extended LB equation, Eq. 
(^), in two and three dimensions and from these compute Hks- We begin with the two 
dimensional case. Here the radius of curvature is a simple scalar, and Eq. ( p^D becomes 

dF ^ dF dF _ 
dt df dp 



rV2 

nav I d(f) cos ( 

l-Tv/2 



(60) 



^ dp' 5 j^p - ^"Tcos,/. j ^(^' ^ ' P' ^) - ^(^' ^' P' 

To find an equilibrium solution, we look for solutions that do not depend upon time, velocity 
direction, or position, and which become the known equilibrium solution for the LB equation 
when the integration over p is carried out. That is we look for solutions F of the form 

F{v,p) = ip,{v)^{p) (61) 

where '{'oiv) = {2nvV)^^6{v — vq) is the normalized, equilibrium spatial and velocity distri- 
bution function for the moving particle with constant speed, which we here denote as vq, 
and confined to a volume V. Then we require that ip{p) be normalized as 
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rfp^(p) = 1. (62) 
It is an easy matter to obtain an equation for ip{p) which reads 

''l^ " -2^0^^ + nav d<P ^ dp' 5 - ^Yacot.^ j ^(pO • (63) 

An inspection of the deha function shows that it vanishes unless p < a/2. Therefore, for 
P > ct/2, we have the simple result 

^(p) = Ae"7 for p>a/2 (64) 

where i = {2nav)~^ is the mean free path length for the moving particle at low density of 
scatterers, and A is a constant to be determined. To treat the distribution function ip{p) for 
smaller values of p we note that we can require that iplp) a.s p 0, since the dynamics 
will increase the value of p during free particle motion and will decrease it to some value, 
still greater than 0, at the instant of a collision. Further we can require that p be continuous 
at p = a/2 since the extended LB equation does not have an explicit delta function of the 
form 6{p — a/2) on the right hand side. Finally we note that the dominant contribution to 
the p' integral on the right hand side of Eq. (|63| ) comes from p' ~ i. Therefore the delta 
function on the right hand side of Eq. ( |63D can be approximated hyf\ 



a cos ( 



^'^(p-(«cos0)/2). (65) 



After inserting this expression in Eq. (^) we find that for p < a/2 



V—- = 2nav 
op 



2(J fOO 



(66) 



where a = ap/2. Since ip is normalized to one, we find thatQ 



more detailed examination of this integral keeping the full delta function shows that the terms 
neglected here are of order nlnn compared to the terms retained, where h = no?. 
^Comparing Eqs. ( p^ ) and (^) one sees that apparently has a discontinuity at p = a/2. 

This, however, is an artifact of the low density approximations we have made. Note that the jump 
in -0 is of relative order no? indeed. 
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( 




2\ V2' 



^(p) = (!/£) 1 



for p < a/2, 



(67) 



and, using the normalization condition on ?/^, we find 



^(p) = (l/£)e"7 for p > a/2. 



(68) 



Combining this expression for ip with Eqs. (^ ^T]) we obtain 



A(+) = h 



•Ks = 2naf (1 — ln2 — C — Inn) for n <^ 1, 



(69) 



in agreement with the result given by Eq. (|26|). 



Now we turn to the three dimensional case. This is somewhat more complicated than the 
two dimensional case since p is a 2 x 2 matrix and not a scalar. However, we can still simplify 
the delta function in the restituting collision integral by noticing that the diagonal elements 
of the curvature matrix grow with time during the free steaming intervals between collisions. 
Consequently, the diagonal elements of p' appearing on the left hand side of Eq. (^6D are 
of the order of the mean free path length immediately before a collision with a scatterer. 
An elementary consideration of the properties of the inverses of 2 x 2 matrices with large 
diagonal elements shows that the dominant contribution to the radius of curvature matrix 
p comes from setting the left hand side of Eq. (|56|) equal to zero. This greatly simplifies the 
delta function appearing in the collision integral on the right hand side of Eq. (|59|). The 
effect of this simplification is that we have neglected terms of relative order a/£, which are 
density corrections to the terms we keep. 

The equilibrium distribution function F{v, p) can be factorized as 



where ^{v) = {Att'^VqV) ^6{v — vq) is the normalized equilibrium distribution function for 
the moving particle. The extended LB equation reduces to 



+na'^v Jq' d(j) Jq"" da J dp[^ J dp'^^ J dp'^i J dp'22 sin cos n»j ^{pij - Pij{4>, a)) ^(p')- 



F{v,p) = ^{v)ij{p), 



(70) 




(71) 
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The matrix elements of p{(f), at) can be obtained by solving Eq. (|56|) under the approxima- 
tion [p']^^ = 0, which can be justified again as a low density approximation by the same 
arguments as in the two dimensional case. The results are 

/ / \ ^ COS ^ / ^ 1 I • 1 \ 

pii{(p,a) = — - — (1 + tan 0sm a) (72 j 

asin2 0sin2a 

Pi2 (0, a) = p2i (0, a) = (73) 

4cos0 

/ / \ d COS (p . 2 / 2 \ /i^a\ 

P22(0,a) = — ^ — (1 + tan 0cos a). (74) 

The restituting term in Eq. ( |7TD contains an integration of ipip') over its arguments, and 
we can set 

^dp'i,{p') = l (75) 

We use if) to compute the average value of Tr(p)^-'^ = [Trp][det p]^-*^, which determines 
the sum of the positive Lyapunov exponents. The pij (0, a) occurring in the delta function 
on the right hand side of Eq. (|7T| ) can be identified at low density with the values of pij right 
after a collision with collision parameters and a. From Eq. ( [73D it follows that these values 
are always equal for pi2 and p2i- Furthermore these quantities do not change in between 
collisions, so we may set ip{p) = /(pii, P225 Pi2)^(pi2 — P2i)- Next we change variables from 
Pii) P225 P12 to pi, P2, P12 where pi, p2 are the eigenvalues of the 2x2 matrix p. The Jacobian 
for the transformation of variables, given by 

J ( P1,P2 \ ^ |Pl -P2I .„f.^ 

[pn,P22j [(pi-P2)2-4pyV2' ^""^ 



will be included in the integrand in Eq. (|80D . After some straightforward algebra we obtain 
the following equation for fi'(pi, P2, P12) = /(pii, P22, P12): 

^ + ^(^1' P12) + '^9{pup2, P12) = 

^ lo^' /o " da sin cos 0| cos 2a\6 (pi - 6 (p2 - (a cos 0)/2) (77) 

^ (Pi2 + (ct COS tan^ sin 2q;) /4) . 
where u = nira'^v is the average collision frequency for the moving particle. The a integration 
can be carried out and we see that g can be written in the form 
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9{pi,p2, P12) = e f 1 - .^k^ilL j h{pi,p2) (78) 

V \p1-p2\J 



where h satisfies 

^ lo^' d<l>^-£^6ip2 - (acos0)/2)5(pi - a(2cos0)-i). 
We can now express the sum of the positive Lyapunov exponents in terms of h as 



(79) 



J dp, J dp2 J dp,, {± + ^) { ,J;:y7t.^,.r^ ) e (1 - Mpi, p.) 

The P12 integration can be carried out easily yielding 



(80) 



max ' ^min 



^ J dpi J dp2 (^y^ + \p2- pl\h{puP2). (81) 



(82) 



We define a new function p{pi, P2) = (7r/2)|p2 — Pi|^(pi) P2) which satisfies 
^ (a?r + &) P2) + MPUP2) = 

^So'^dct^'^ - COS0) 5{p2 - (acos0)/2)5(pi - a(2cos0)-i). 

Introducing pi(p) = dp2p{p, P2) and P2{p) = /o°° dpip{pi, p) enables us to express the sum 
of the Lyapunov exponents very simply as 

roo \ 

Knax + K^^n =11 ^p- Pi{p) (83) 

where the Pi{p) satisfy 

dpi(p) /"^/^ 
V — \-upi{p)=2vj d(j) sin (j) cos (f)6 {p — {a cos (f))/ 2) and (84) 

3p2{p) /"""/^ 

V — Vijp2{p)=2ul d(f) sin (p COS (l)5{p — a(2 COS (p)^^). (85) 

CJ f) J 



Solving Eqs. (|8^ , ^) for ]9i,p2 and inserting the solution into Eq. (|83| ) one obtains an 
expression for the sum of the positive Lyapunov exponents which agrees with that of Eq. 
(p8|). It is worth mentioning that we can solve Eq. (ff^) to provide g{pi, P2, P12) as an 
explicit function of the variables pi,P2,Pi2 using the method of characteristics. As this 
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solution will prove useful in subsequent papers, we outline the method in Appendix B. 
There we also briefly indicate how the individual Lyapunov exponents can be calculated 
using g{pi, P2, P12) and simple results from the theory for eigenvalues of products of random 
matrices ||23|. While the results obtained with the extended LB equation properly agree with 



those obtained by more direct kinetic theory arguments, we will need to use the extended 
LB equation in order to obtain the Lyapunov exponents and KS entropies for the spatially 
inhomogeneous systems that occur when one considers escape-rate methods for connecting 
chaotic quantities with transport coefficients. This will be treated in a subsequent paper. 

VI. THE NEGATIVE LYAPUNOV EXPONENTS AND THE ANTI-LORENTZ 

BOLTZMANN EQUATION. 

We now turn our attention to the Lyapunov exponents that characterize the exponential 
convergence of trajectories on a stable manifold in the 2d — 1 dimensional constant energy 
surface in the phase space of the moving particle. We recall that two arbitrary but inflnitesi- 
mally nearby trajectories will certainly separate eventually with time. We have used this fact 
to derive formulas and explicit expressions for the positive Lyapunov exponents. However 
Liouville's theorem showing that the measure of a small region of phase space is constant in 
time as one follows the motion of points initially in that region implies that there must be 
a compensating set of negative Lyapunov exponents which act in concert with the positive 
ones to keep phase space measures constant in time. Moreover, the fact that the Lorentz gas 
is a symplectic Hamiltonian system has as a consequence the existence of a conjugate pairing 



rule [0 . That is for such a system, the Lyapunov exponents come in positive and negative 
pairs such that the sum of each corresponding pair is zero. Thus, in this case at least, the 
calculation of the negative Lyapunov exponents is trivial - they are just the opposites of the 
positive ones. However, in a subsequent paper where we treat thermostatted systems, this 
form of the conjugate pairing rule no longer holds p8|-|30|] and we will need to flnd methods 



to compute both positive and negative exponents individually. 
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The most obvious way to obtain the negative Lyapunov exponents is to compute the 
positive Lyapunov exponents for the time reversed motion. Upon time reversal trajectories 
that approach each other in the forward motion will separate. In fact almost all trajectories 
will separate in both the forward and the time reversed direction, but in general (i.e. for non- 
symplectic systems) with different exponents. In the forward motion they will separate with 
rates given by the positive Lyapunov exponents, and in the time reversed motion with rates 
equal to the magnitudes of the negative Lyapunov exponents. Thus to calculate the negative 
exponents we consider the binary collision dynamics already discussed before, but look at the 
time reversed motion. If, in the forward time direction, the moving particle is uncorrelated 
with a scatterer before collision, in the time reversed motion it will be uncorrelated with the 
same scatterer after the coUision. We therefore should consider a kind of backwards kinetic 
theory, where the particles are uncorrelated with the scatterers after their collisions instead 
of before them. This will differ in important respects from the ordinary Lorentz-Boltzmann 
equation. In order to illustrate this we consider first a calculation of the sum of the negative 
Lyapunov exponents for the 2-d and 3-d dilute Lorentz gases, using a simple kinetic theory 
argument similar to that used in Section II. 

A. A Simple Kinetic Theory Method for the Sum of the Negative Lyapunov 

Exponents 

We first consider the 2-d case. We wish to follow an infinitesimal trajectory bundle which 
is contracting and remains contracting for all future times. Such a bundle is illustrated in 
Fig. 1. We follow the motion of this bundle from scatterer 1 to scatterer 2. In order that the 
bundle remain contracting after the collision with scatterer 2 we require that the radius of 
curvature of this bundle be very close (with corrections of order na^) to the value a cos 0/2 
just before the collision with scatterer 2. We denote the direction of the velocity of the 
moving particle just after the collision with scatterer 1 by the angle 9 with respect to some 
space fixed axis, and compute the negative Lyapunov exponent in the following way. Since 
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we know that the radius of curvature just before the colhsion with 2 is a cos 0/2, the radius 
of curvature just after the colhsion with 1 must be a cos(j)/2 + vt, where t is the time interval 
between the collision with 1 and 2. From this it follows that 



V f dT[p{T)]-^ = In 
Jo 



2vt 



acos( 



(86) 



where p(r) is the radius of curvature of the contracting bundle at times < r < t before 
the collision with 2, and vt ^ a. The time average of this expression corresponds to the 
result obtained by combining Eqs. ([T7|) and (|2l|) , in agreement with the conjugate pairing 



rule for symplectic systems. In the case of a thermostatted system however, Eqs. (pif ) 
and (|8^) are replaced by expressions that depend on the velocity angle 9 in different ways. 
As a consequence the positive and negative Lyapunov exponent are no longer each others 
opposite. 

The three dimensional case proceeds in exactly the same way. We follow an infinitesimal 
contracting trajectory bundle from a collision with scatterer 1 to a collision with scatterer 
2, such that before the collision with scatterer 2, the radius of curvature matrix is given by 
[(2 cos0P)/a]^^ where P is defined in Eq. (^T]). We then find easily that if t is the time 
interval between the collisions of the moving particle with scatterers 1 and 2, then the radius 
of curvature matrix at some time r between zero and t after the collision with scatterer 1 is 

For calculating the sum of the positive Lyapunov exponents one needs the average again 
of Tr[p(r)]~^, but now under the Stosszahlansatz for the postcollisional (under the time 
reversed dynamics, so in reality the precollisional) coordinates. However, due to the time 
reversal symmetry of the dynamics the distribution of the postcollisional coordinates (f)' and 
a' is the same as that of and a, and of course the distribution of intercollisional times also 
is the same for forward and backward motion. Therefore the averaging procedure yields the 
same results as for Eq. (^Op and the sum of the negative Lyapunov exponents is given by 
the opposite of Eq. (p8|). 
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The individual Lyapunov exponents can be obtained using results from the theory for 
eigenvalues of products of random matrices, as described in Appendix B. The results are, 
as expected, that each negative exponent is paired with a positive one such that their sum 
is zero. Now we turn to the method of distribution functions for calculating the negative 
Lyapunov exponents. 

B. The Extended Anti-Lorentz-Boltzmann Equation 

In order to use distribution functions to compute the negative Lyapunov exponents, or 
their sum, we need to construct a Boltzmann-like equation for the time reversed motion. If 
one reviews the derivation of the Boltzmann equation, one sees that the colliding particles 
are taken to be uncorrelated before the direct or the restituting collisions. If we were to 
look at the time reversed motion, the terms "before" and "after" are interchanged, and for 
the time reversed motion the colliding particles are uncorrelated after the collisions rather 
than before. That is, referring to Fig. 1 again, the time reversed motion has a collision of 
the moving particle with scatterer 2, followed by a collision of the particle with scatterer 
1. Before the collision with 2, the moving particle is correlated with 2, but not so after 
the collision. After the collision, the radius of curvature or the radius of curvature matrix 
is given by the values acos0/2 in two dimensions, or [2cos0P/a]~^ in three dimensions. 
The radius of curvature (2d) or the diagonal elements of the matrix (3d) grow over the time 
interval t until the collision with 1. Since the moving particle is correlated with the scatterers 
before the collisions and not after them, when deriving a Boltzmann-like equation for the 
distribution function, we must apply the Stosszahlansatz to the exiting collision cylinders 
after the collisions and not to the entering cyhnders before coUision, as is the usual case. 

In order to clarify this procedure we first consider a derivation of the anti-Lorentz- 
Boltzmann equation (ALBE) for the distribution function /_ (r , v, t) for the position and 
velocity of the moving particle at time t. We will then generalize this derivation by adding 
the radius of curvature variables. Since there is no particular problem with the time reversal 
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of the free streaming of the particle, we can write the ALBE in the form 

^^#^ + -^ = r^-r- (88) 

ot or 

where Fl represents the rate at which particles are produced at r with velocity v and Fl 
represents the rate at which such particles are lost. Noting the fact that we should apply 
the Stosszahlansatz to the exiting collision cylinders, we see that the rate at which particles 
with velocity v are produced is 

Tt5r5m = na'^-^ J dn{h ■ v)f-{f, v, t)SrSvSt (89) 

where n ■ v > 0. That is, particles with velocity v' = v — 2(n ■ w)?! collide with scatterers 
and produce particles with velocity v, but the distribution function we need to calculate the 
rate of these collisions for is not /_(r, -y', t), but /_(r, -u, t). Similarly 

r: = na'^-^ J dn{n ■ v)f-{f, v, t) (90) 

where {"h ■ v') > 0, since particles with velocity v' are produced when a particle of velocity 
V collides with a scatterer with collision vector n. Putting these terms together, we obtain 

^I^^^^ + v-^ = na'-' Jdh\v. h\ [/_(f, V, t) - /_(f, v\ t)] (91) 

where we have used the fact that |t?-ri| = \v' ■ h\, and the integration in Eq. (^Tj) is over 
a semicircle (for (i = 2) or a hemisphere (for = 3). This equation looks exactly like 
the Lorentz-Boltzmann equation, Eq. (^7]) , except for the fact that the collision integral 
has the opposite signs in Eqs. (|9T|) and (|^) |31[|. The ALBE, Eq. (|9TD has the usual 



equilibrium distribution function as a stationary solution, although it is a highly unstable 
solution, any deviation will tend to grow exponentially in time. In fact, for this reason it is 
an ill-posed equation, as arbitrarily small initial deviations can grow at arbitrarily large rates 
and therefore its solution is not well-defined. This ambiguity can be removed by requiring 
that the actual solution describing a physical system after a long enough time that all these 
rapidly decaying solutions (observed in the forward time direction) have died out, is the time 
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reverse of the solution of the ordinary LBE (in a closed isolated system or system coupled 
to a single heat bath this will become the equilibrium solution for long times). So there is 
little use in employing the ALBE as such, since to obtain its physically relevant solution one 
has to solve the ordinary LBE anyway. The extension of the ALBE discussed in the next 
paragraphs on the other hand does provide a useful tool for calculating negative Lyapunov 
exponents. 

To this end we need to include again the radius of curvature matrix as variables in the 
distribution function, F_{f,v,p,t). The equation for F_ is constructed as before, i. e., we 
write 



rl(p)-r:(p). 



(92) 



dt dr ^ dpi. 

To compute r+(/3) we multiply the rate at which particles with velocity v are produced no 
matter what their radius of curvature matrix might be, with the fraction of particles with 
velocity which produce particles with radius of curvature p after collision. That is 



(93) 



rl (p) = na"^ ^ / dh\v ■ h\ J dp'F^{f, v, p', t) x 
/ dp"F_ (f, V', p", t)5ip - Pip"))] [j dp"F_ (f, V', p", t)] . 

Similarly, to compute ri(p) we multiply the rate at which particles with velocity v' are 
produced due to collisions of particles of velocity v with scatterers by the fraction of particles 
of velocity v that have the radius of curvature matrix p. That is 



(94) 



r:(p) = na'^-^ Jdn\n-v\Jdp'F4f,v',p',t)F(f,v,p,t)x 

[Jdp'F_{f,v,p',t)]-\ 

We can now assemble all of these results into an extended anti-Lorentz-Boltzmann equation 
(EALBE) which reads 

dF^ ^ dF_ t^dF_ 



dt dr fr{ dpii 

^a'^~^ j dh\v-h\ j dp'F_{f,v,p',t) x 

J dp"F_{f,v',p",t)6{p-p{p"))'\ [I dp"F_{f,v',p",t) 
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J dn\n ■ v\ J dp'F_{f,v', p',t)F{f,v, p,t) x 



J dp'F.{f,v,p',t) 



(95) 



Here, too, n is integrated over a semicircle or a hemisphere. This comphcated and non- 
hnear looking equation (in fact one should first solve the ALBE, substitute its solution for 
/ dp'F^{f, v', p', t) in Eq. (^) and then solve the EALBE; both equations to be solved then 
are linear) will form the basis of the calculation of negative Lyapunov exponents and their 
sums in the more complicated cases to be studied in subsequent papers. We should note here 
that without further conditions this equation, like the ALBE discussed above, constitutes an 
ill-posed problem. The additional condition that regularizes its solution is that the integral 
over p of the distribution function yields the time reverse of the solution of the ordinary 
LBE. In the equilibrium case considered here the equation simplifies enormously. In fact if 
we use the condition that the system is in a spatially homogeneous equilibrium state, that 
the integral of the distribution function F_ over all elements of the radius of curvature matrix 
is the equilibrium distribution function which is independent of the velocity direction, 
we immediately obtain Eq. (7T) for the spatially homogeneous case. Therefore for this 



equilibrium situation, F_ produces the same Lyapunov exponents as in the forward motion, 
with the exception of the appropriate change of sign due to the time reversed nature of the 
motion considered here. 

VII. COMPARISON WITH SIMULATIONS AND DISCUSSION 

We have found that the quantities calculated in this paper, the Lyapunov exponents and 
the KS entropies, expressed as functions of the dimensionless density n, have the general 
form 

K = Au[-\nh + B + o{l)], (96) 
hKs = Au[-\nh + B + o{l)], (97) 
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where v and h are the coUision frequency and reduced density, given by z/ = 2nav,h = na^ 
in two dimensions and v = irna'^v, h = nna^ in three dimensions, and A, B are constants 
which we have determined. In Table I we compare the theoretical results for A and B 



with values for the same coefficients, as obtained by Dellago and Posch from computer 
simulations of two and three dimensional hard sphere Lorentz gases. The results are in 
excellent agreement for the coefficient A and there are minor discrepancies in the B values, 
probably due to the fact that the simulation analysis is difficult at the low densities where 
the theoretical analysis given here applies. 

It is remarkable that the Lyapunov exponents for the 3 dimensional completely isotropic 
random Lorentz gas are different at all. In fact, on the basis of the results obtained in 
leading order in the density, it has been conjectured that all positive as well as all 
negative Lyapunov exponents are equal. The methods used in our approach allow a very 
transparent explanation for the differences, now also confirmed numerically, between the 
Lyapunov exponents. Especially it becomes clear why all Lyapunov exponents coincide 
in leading order but differ in next to leading order. The reason for this is related to the 



different nature of the terms contributing to the different orders |]TE[ . The terms proportional 
to nlog(n/2) + C in Eqs. ( |27| , result from averaging over functions, which only depend 
on the time of free flight. Due to the isotropy of the free flight this has to lead to equal 
Lyapunov exponents. The differences arise on averaging over functions, which depend on the 
collision parameters and a. To understand how the scattering process, which is isotropic 
for a single trajectory, can cause the Lyapunov exponents to be different, it is worthwhile 
considering scattering events in more detail. 

For calculating the Lyapunov exponents we have to analyze the scattering of two close 
by trajectories. Therefore it becomes possible for a given scattering angle (j) to distinguish 
whether 5rx is in the plane V spanned by the normal h on the sphere and the impact 
velocity v or perpendicular to V. In this way the isotropy of the scattering process becomes 
effectively broken. This is reflected in the eigenvalue structure of the radius of curvature 
matrix = {2 cos (f)P / a)~^ (see Eq. (|57D). There are two eigenvalues pi = a cos 0/2 and 
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P2 = a/(2cos0). Here the eigenvalue pi corresponds to the eigendirection ei which is in 
the plane "P, and the eigenvalue p2 corresponds to the eigendirection which is perpendicular 
to the plane P, with both eigendirections perpendicular to v. Note that for = the 
eigenvalues are the same. This can be understood by realizing that for = 0, the particle 
hits the sphere head on , i.e v is parallel to n. For this special case the two eigendirections 
are clearly equivalent, but for other values of 0, this symmetry is lost, and the eigenvalues 
differ. Therefore we may conclude that the lack of degeneracy of the positive Lyapunov 
exponents is due to the lack of rotational symmetry when the nearby trajectories hit the 
sphere. 

We conclude with a number of remarks: 

1) The results given here can be extended to higher densities in a number of ways. In 
particular, BBGKY hierarchy methods are being developed to provide a systematic den- 
sity expansion of the Lyapunov exponents and the KS entropies beyond the low density 
results obtained here. This will be especially important when non equilibrium situations are 
considered, since there one may see the effects of long time tail phenomena on the chaotic 
properties of the system. 

2) As remarked earlier, we relied upon the low dimensionality of our systems to obtain 
all of the relevant Lyapunov exponents. For a four dimensional Lorentz gas we would need 
to use more sophisticated techniques to obtain all of the Lyapunov exponents, since the 
methods given here could only provide values for the largest exponent and the sum of all 
of three of the positive exponents. We could not then resolve the two smaller positive 
exponents, but only could get their sum. 

3) We have not analyzed here a particularly interesting quantity that gives a more general 
characterization of the chaotic properties of the system, namely the Ruelle, or topological, 
pressure • This quantity has the formal structure of an equilibrium free energy and it 
depends upon a temperature-like parameter, (3. The results obtained here characterize the 
chaotic properties in the neighborhood oi (3 = 1. Elsewhere it has been shown that for the 
Lorentz gas, where the disorder is static, in the thermodynamic limit the topological pressure 
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exhibits a localization transition as a function of (3, it is dominated by contributions from 
particles which are localized within the largest dense cluster for /5 < 1, and in the largest 



region without scatterers for [3 > \ |3^. It would be very valuable to obtain these results 
using kinetic theory methods in addition to the more rigorous analysis given in Ref. Q . 

4) It should be straightforward to extend the results obtained here to other potentials of 
interaction between the moving particle and the fixed scatterers, at least for dilute systems. 
It is well known that the Boltzmann transport equation can be applied to gases that interact 
with other than hard core potentials Jl^. In fact a wide variety of interaction potentials 
may be used in the Boltzmann equation to determine the transport properties of the corre- 
sponding gases. Certainly Lyapunov exponents can be calculated for Lorentz gases where 
the moving particle has other than hard core interactions with the scatterers, as well. 

5) Finally we mention that the method given here can be adapted to cases where all the 
particles in the system are moving, namely, a dilute gas. Expressions for the KS entropy 
p5| and the largest Lyapunov exponent of two and three dimensional gases with short range 
forces have already been obtained this way. 

In the next paper we will extend the results here to open systems with escape and 
compute the escape rates, Lyapunov exponents, and KS entropies characterizing the fractal 
repeller that underlies diffusion in open systems with absorbing boundaries. This work 
will make heavy use of the extended LB equation, and we regard the present paper as an 
introduction to the next one. 
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TABLES 



QUANTITY 


THEORY 


SIMUL. 


2d: A+ 






A 


1 


0.995 ± 0.009 


B 


0.423 


0.463 ± 0.083 


OA. \ + 






A 


1 


0.990 ± 0.089 


B 


0.309 


0.387 ± 0.7-lG 








A 


1 


0.992 ± 0.084 


B 


-0.077 


-0.015 ± 0.715 


3d: hKS 






A 


2 


1.982 ± 0.173 


B 


0.166 


0.372 ± 1.461 



TABLE I. Comparison of Theoretical and Simulation Results 
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APPENDIX A: BINARY COLLISION DYNAMICS AND THE RADIUS OF 

CURVATURE MATRIX 



Here we give a brief review of the derivation of the formulae for the change in the spatial 
and velocity deviations at a binary collision used in Eqs. (|[ 0) based on the method of 
Dellago, Posch and Hoover We relate these formulas to the expression used for the 

change in the radius of curvature matrix, given by Eq. (|35D , which, in turn, was discussed 



by Gaspard and Dorfman Consider a trajectory of the particle moving among the 

scatterers. We denote the initial position and velocity of this trajectory by a;(0), and the 
position and velocity at time t later by x{t). Consider also a trajectory that is obtained 
by an infinitesimal displacement of the initial position and velocity to x(0) + 6x{0), and 
denote the position and velocity of the second trajectory at time t by x(t) + 6x(t). We 
require the two trajectories to be infinitesimally close. Our goal is deriving equations for 
6x{t) = {6f{t),6v{t)). We take all trajectories with the same energy, which leads to the 
condition that 6v{t) ■ v{t) = 0, so the velocity deviation is always perpendicular to the 
velocity v(t). We can also set the position deviation 6r{t) ■ v{t) = 0, since this simply 
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requires that the position deviation is perpendicular to the velocity at the initial time. 

Now in between collisions 6f{t),6v{t) satisfy Eqs. (|^, However the change in these 
quantities at collision is more complicated. To analyze this change we suppose that the 
trajectory with x{t) has a collision with some scatterer at time r. Then immediately after 
the collision the velocity has changed to = v — 2{v ■ n)h where n is a unit vector in 
the direction from the center of the scatterer to the point of contact, and v is the velocity 
immediately before collision. The displaced trajectory will have a collision at a slightly 
displaced time r + (5r, and at a slightly different point on the same scatterer located by unit 
vector h + with n ■ 5n = 0. By examining the scattering equations for the displaced 
trajectory one easily finds 

5v^ = (1 — 2hh) ■ 6v — 2[{v ■ n)6n + {v ■ 6'h)n] (Al) 

5f(r) + 6tv = a6h (A2) 
Sr^ = -Stv^ + a5h. (A3) 

We use the condition n ■ 6h = to obtain 6t = —{n ■ 6f{T))/{v ■ n). Simple algebra leads to 
Eqs. in the text. It is important to note that 5r(r) is the spatial deviation when the 

"main" trajectory has a collision, while 5r^ is the spatial separation at the instant that the 
displaced trajectory has a collision. (It is easy to visualize this in the case that 5t > 0, as 
illustrated in Fig. 2.) 

Now we introduce the radius of curvature matrix as given by Eq. (|33|) . We note that p 
is a matrix of rank d — 1. If we substitute the definition of this matrix, Eq. (|33|) into Eq. 
(I), we obtain 

p+-Sv^ = V- p- ■ Sv-, (A4) 

where U = U^^ = 1 — 2nn is a reflection matrix and is discussed in Ref. [|19|. Here 
the superscripts +, — denote after and before collision respectively. The operators U have 
determinant -1 and they ensure the proper orientations, at collision, of the planes in which 
the radius of curvature matrices are defined. We now substitute Eq. (^ for 6v'^ into Eq. 
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use the relation Eq. (^), and some elementary matrix manipulations to obtain Eq. 
(p5|). As mentioned in the main text the inverse matrices [p^^^*-^-*] and [p^^'-^-*] are defined 
in the subspaces orthogonal to v and v' respectively. 

We can obtain a considerable simplification of the analysis of the spatial deviation vector, 
in Eq. by using the properties of the matrices U. That is, we can easily show that 
it is possible to express the spatial deviation vector in terms of radius of curvature matrices 
p all of which are defined in a plane perpendicular to the initial velocity ^^(0). To see this, 
consider the right hand side of Eq. (|37|) for the case that n = 2. We have 



(A5) 

[l±(0)+t;ri,op-i(0)]-5fx(0). 

Here l±{i) is a unit operator in the plane perpendicular to the velocity vector after the i — th 
collision. Note that Eq. (|35|) allows us to write 

Pu'^^\h) = V{l)--pZ'^^\h).V{l), (A6) 

where 

. 9 [ 1 1 

(A7) 



(Jj 



v{0)ni + niv{0) - ^ , . fiihi - {v{0) ■ hi)l 
(f(O)-ni) 



It is important to note that the operator defined by the square bracket appearing on the 
right hand side of Eq. { \K^ ) is orthogonal from the right and left to "^(0), as is p^. Therefore 
the operator p^^^ is to be evaluated in the plane perpendicular to "^(0). Furthermore the 
unit operator 1±(1) appearing in the same bracket with Pu^^^\ti) can be written as U(l) ■ 
1^(0) • U(l). From this it follows that 

6/-\t,) = U(l) ■ [U(0) + vr,yp-'^^\t,)] . [1^(0) + p-\0)] ■ 5r-l(0). (A8) 

If this procedure is followed through each successive collision, Eq. (|37| ) can easily be written 

as 

5ri-)(t„) = V{n - 1) • V{n - 2) ■ ■ -Vil) ■ 5r=i"^(t„), (A9) 
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where 

6^-\Q = [U(0) + p-'^^\t,.-i)] ■ ■ ■ [U(0) + p„-^(0)] ■ Sf^iO). (AlO) 

The product of the U matrices appearing on the right hand side of Eq. (A9) have determi- 
nant ±1, of course, and have no bearing on the exponential growth of the spatial deviation 
vector. As a result one can carry out all calculations of the Lyapunov exponents and the 
KS entropy in a coordinate system defined in the plane perpendicular to the initial velocity 
■^(0). As a result, all of the U operators can be dispensed with in the calculation of the 
Lyapunov exponents provided one uses operators as well as unit operators 1±(0). This 
procedure "unwinds" the trajectory of the moving particle and allows all collisions with the 
scatterers to be treated in one coordinate system. This is useful when one wants to avoid 
neglecting the left hand side of Eq. (^61), to treat systems at higher densities, and/or to use 
the methods of random matrix theory in a convenient way, as we do in the next appendix. 

APPENDIX B: THE SOLUTION OF THE EXTENDED LB EQUATION USING 
THE METHOD OF CHARACTERISTICS, AND RANDOM MATRIX METHODS 

Here we indicate how Eq. ([7ll) can be solved as a differential equation in three variables 
using the method of characteristics ||3^. We begin with Eqs. ([77D, and change variables to 



(Ji = 2pi/a, (T2 = 2p2/a,ai2 = 2pi2/a. We also set g{pi, p2, Pu) = 8/a^G{ai, a2, ctu) and 
introduce the scaled density n = nira^. Then the equation for G is 



4i^e(i - a,)e (1 - ^) {^)s{a2-j-^). 



where K = To use the method of characteristics we write 



\dai da2 J ds 

where 



« +JL)g = #G (B2) 
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£ai = 1; ai=a^ + s 

= 1; a2 = aO + s . (B3) 
Aau = 0; au = cr^u- 

This substitution converts the partial differential equation into a simple differential equation 
which can be solved by elementary means as a function of s, once we have specified the 
appropriate boundary conditions. Further, in the cti, o"2 plane, the s integration corresponds 
to an integration along a set of lines given by Eqs. (P3|), called the characteristic lines. The 
inhomogeneous term on the right hand side of Eq. ( [Bl| ) is zero everywhere except on the 
line ^2 = l/a"i; (Ti < 1. We look for solutions that vanish at Ui = and at (T2 = 0, and we 
note that the solution G of the homogeneous equation, expressed in terms of s, has the form 

G = Goe-^. (B4) 

All of the conditions can be satisfied if G vanishes in the ai,o'2 plane, except in the region 
defined by the curves cr2 > (Xi, ai > 1 and (J2 > 1/o-i,(Ti < 1. See Fig. 3. Choosing 
cr° = z < 1, cr2 = V^i ^^'^ ^12 = ^^12 •= sin (2d) |a2 — ai|/2 the characteristic lines are given 
by 

ai = z + s; z < 1 

CT2 = l + S. (B5) 

(712 = sin(2d) 

We then arrive at the equation 

+ = K&{1 - z)&{27r - a)&{&) (^-^^ (^^j 6{s). (B6) 

The probability density for the random variables s,z,a is given by f{s,z,a) = 
J{s, z, a)G{s, z,a) where J is the Jacobian J = |(9(cri, (72, cri2)/<9(s, 2;, d)| = (1 — z^){l + 
z^)/z'^. f obeys the equation 

^J + ^f = KQ{1 - z)e{2n - a)Q{a)z6{s). (B7) 
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This equation can now be solved along the characteristic lines in the region indicated in Fig. 
3. It has the simple solution 

f = —9(1 - 2)e(s)e(a)e(27r - d)zexp(--s). (B8) 

When the solution is inserted in the expression for h^s = v{Tr\/ p), 

f 1 1 2f /■ - 11 

hKs = v / dpidp2dpug{pi, p2, Pi2){ \ ) = — / dsdzdaf{s,z,a){ — ■ h — ^ — ■ — ), 

J Pi P2 a J z + s 1/z + s 

(B9) 

the results obtained agree with those given by Eq. (|28|) . 

The calculation of the maximum Lyapunov exponent starts from the observation that the 
deviations Sf±{t) are given by a product of random matrices acting on the initial deviation 
Sr±{0) (see Eq. (0)). The random matrix R(tj+i,tj), Eq. (^51) , depends on the time of 
free flight between two collisions and the initial matrix p+ given by ((a(2 cos0))P~^. P 
(see. Eq.(^)), is parameterized by two collision angles 0, a. Another way of parameterizing 
these initial matrix is by using its random eigenvalues z := = cos0, = 1/-2 and the the 
random angle variable a. Notice that the characteristic lines are identical to the free flight 
solutions for the eigenvalues of {2/a)p and its off diagonal element ai2 (Compare eq. ([B3|)). 
They are chosen such, that the initial conditions of the characteristic lines {s = 0) are the 
values of <Ji,ai2 immediately after a scattering event with scattering angles (p,a, Therefore 
( |B5D can be interpreted as the distribution function for the random variables z = cos0 
(cosine of the polar scattering angle), a (azimuthal scattering angle) and the time of free 
flight s. 

The calculation of the maximum eigenvalue is now a straight forward application of the 



theory of products of random matrices p3|. It follows from Eq. (|4^) that the maximum 
Lyapunov exponent is given by 

lim -log|nRMO)|/|MO)| (BIO) 

t—>oo t 

with IIR = rijlo-f^O'i + l)- Using an analogous decomposition of the product as in Eq. 
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(p!T|) and the identity t = J2iLi ■^j? valid right after the Nth colhsion, with si the time of free 
flight between collisions i — 1 and i, ( P10|) is equivalent to 



Xmax = j^{^og\R{z,a,s) ■ e^\) . (Bll) 



Here is a unit vector = (cos V', simp) and (. . .) indicates an average over the distribu- 
tion of {z,a,s) plus an additional average over a stationary distribution D{iIj) of directions 
e^. This distribution is a solution of a Frobenius-Perron equation. 

D{^|J) = r d^'Di^')5i^' - cos-^( '"?'^^^'^'^;'^:f^ ) (B12) 
JO |-r{,(z, a, s) ■ e^l 

As we will show shortly this additional average is not necessary in the equilibrium case since 
ip can be absorbed in a redefinition of the azimuthal scattering angle a. The equation for 
R can be derived from the equation (|45|). 

^ = VP (B13) 
as 

In the equilibrium this is easily solved using p(s) = p"*" + v 1 s 

R(^, a,s) = l + p+-\s ^ s. (B14) 

a 

The last approximation is valid in the low density limit, since the most important contribu- 
tions to the average over s come from large time of free flight s. P is defined in Eq. 
With this we obtain 



\R{z, a, s) ■ eVI = \ {—^Y + ( + 2 ^ ^ ^ cos(2(a - (B15) 

V z z z-^ 

From Eq. (|B15| ) it is obvious that |R(z, ■ e^| is statistically independent of il) i.e. the 
dependence on -i/; can be absorbed in a redefinition of a. To obtain Xmax we have to calculate 
the average of Eq.( |B15| ) using the distribution function Eq. (PH]). The result agrees with 



Eq. 
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FIGURE CAPTIONS 

Figure 1. A trajectory from scatterer 1 to scatterer 2 with bundles of expanding and 
contracting trajectories indicated. 

Figure 2. The arrangement of the position and velocity vectors for one trajectory at 
the instant of a collision of the moving particle with the scatterer and at the instant of the 
collision of the particle with the same scatterer along an infinitesimally displaced trajectory. 

Figure 3. The solution regions for Eq. ( |B1| ) using the method of characteristics. A 
characteristic line is indicated parallel to the line di = cr2- 
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